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Abstract 

This work gives an algorithm for computing the degrees of freedom of estimators 
of Allan and Hadamard variances, including their modified versions. A consistent 
approach is used throughout. 


1 INTRODUCTION 

This work gives an algorithm by which one can compute error bars for measurements of frequency 
stability variances in the presence of power-law phase noises. These stability variances fall into two 
categories: unmodified variances, which use dth differences of phase samples for d = 2 (Allan) or 
3 (Hadamard), and modified variances, which use dth differences of averaged phase samples. The 
corresponding sampling functions that act on phase and frequency are shown in [1], Each variance 
is defined as a scaling factor times the expected value of the squared differences. Unbiased estimates 
of this variance are computed from available phase data by taking time averages of the squared 
differences. The usual choices for the estimation stride (the time step) are the sample period To 
and the averaging period r, a multiple of to- These give, respectively, the overlapped estimator 
(OE) and nonoverlapped estimator (NOE) of the stability variance (although “nonoverlapped” is 
a misnomer; there is always some overlap between summands). 

The new algorithm, which computes the equivalent degrees of freedom (edf) of a variance estimator, 
can replace several ones currently in use with a single, complete, and consistent method. Specifically, 
this algorithm covers the OE and NOE of the unmodified and modified Allan and Hadamard 
variances for all common noise types (—4 < a < 2) at all applicable sample sizes and averaging 
factors. Previously, for the NOE of Allan and Hadamard variances, 1-sigma confidence limits were 
generally set by scaling the measured deviation by a noise-dependent factor K a /y[M, where M is 
the number of summands ([ 2 - 6 ], For the OE of Allan variance, empirical edf approximations ([ 7 - 
10 ]) were generally used along with chi-squared statistics; non-empirical methods for both Allan 
variance estimators were also published ([ 11 - 14 ]). Although an edf computation for the OE of 
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modified Allan (and time) variance was first given in [ 15 ], the simpler approach of [ 16 ] was gen- 
erally used; an alternative unpublished approach using the Hadamard edf formed the basis of the 
new algorithm. As an example of the application of the new algorithm, the Stable32 program for 
frequency stability analysis [ 17 ] has adopted it (since Version 1.41) instead of the multiple algo- 
rithms previously used for setting confidence limits and error bars. 

We wish to maintain a clear distinction between a stability variance and its estimators; in particular, 
we treat both the OE and the NOE of each variance. Even though the OE usually has lower 
uncertainty than the NOE, the NOE is convenient when phase data are processed in real time or 
read sequentially from a file. Indeed, the TSC 5110A Time Interval Analyzer [ 18 ], in its “averaging” 
mode, computes an NOE of modified Allan variance. 

The goal of the new algorithm is not closed formulas, but fast numerical results whose accuracy is 
adequate for the purpose at hand. All the calculations are based on the same theoretical principles, 
with no empirical formulas. For each r, one must choose a dominant phase-noise power law, 
S x (/) ~ Cf a ~ 2 , where a £ {2, 1, 0, —1, —2, —3, —4} (white PM to random run FM); see [19] for a 
method of power-law identification. The phase noise is assumed to be approximately bandlimited 
to the Nyquist frequency 1/ (2ro). 

Not covered are effects of trend removal, drift removal for Allan variance in particular; the dth phase 
differences are assumed to have zero mean. One can use Hadamard variance to obtain stability 
results that are invariant to linear frequency drift. Special long-term stability statistics, such as 
total Allan variance [20], total Hadamard variance [19], and Theol [21], are not covered; these 
require their own treatments. 


2 THEORY OF OPERATION 


Although the presentation of the algorithm is self-contained, we give a brief account of the theory 
behind it. The algorithm’s output is the equivalent degrees of freedom (edf) of an unbiased estimator 
V of a stability variance a 2 = EV. Define 


v = edf V 


2 (EV) 2 
var V 


2 a 4 
var V 


(1) 


It has been observed empirically (but not exhaustively) for these estimators that (i'/cx 2 ) V has ap- 
proximately a xt distribution. Thus, having computed v and observed V, one can obtain confidence 
intervals of form vV/x^ < a 2 < vVjxy from xt levels x\ < x '2 [ 7 ]. 


The model for phase x (t) is the ro-difference of a pure power-law process: 

x (t) = A To w (t) , (2) 

where w (t) is a continuous-time process with spectral density Cf a ~ 4 for all / > 0, and A is the 
backward difference operator. Then S x (/) is asymptotically proportional to f a ~ 2 as / — > 0 and 
has approximate bandwidth 1/ (2 To); this is the first reason for using w (t). 


Now let 

z (t) = A x A £ w (t ) , 


(3) 
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where e = tq or r. For e = To we have z (t) = A (t), which leads to an unmodified variance. For 
e = t = rriTo we observe from (2) that 

m— 1 m — 1 

A r w ( t ) = w (t) — w (t — m,To) = A To w (t — htq) = x(t — utq) = rnx (t ) , 

n = 0 n = 0 

where x (t) is a discrete-time average of m samples of x. In this case, z[t) = mA^x(t), which 
leads to a modified variance [16]; this is the second reason for using w (f ) . In either case we as- 
sume that z (t) is a stationary zero-mean Gaussian process with autocovariance (ACV) function 
s z ( t ) = E z (u + t) z ( u ). 


Ignoring the conventional scaling factors, we define the stability variance and its estimator by 

1 M 

a 2 = Ez 2 (t) , V = — z 2 ( nS ) 7 ( 4 ) 

n= 1 

where the stride <5 is To for the OE and r for the NOE. The number of terms M depends on the 
estimator type and the number of data. We have EV = a 2 = s z (0). Then cov \z 2 (t) ,z 2 (u)] = 
2s^ (u — t), and 


var V = 


2 

w 


^(( n 2 


ni ,722 = 1 

The definition (1), after substitution of (5), simplifies to 


ni) 5 ) . 


1 

edf V 


1 

M 


1 + 


M— 1 

*W) § 



( 5 ) 


( 6 ) 


The ACV s z (t) is obtained from (3) by applying a difference operator of order 2d + 2 to the 
generalized autocovariance (GACV) s w (t) of the power-law process w (t) [13,16]: 

Sz (t) = (A r A_ T ) d A £ A- £ s w (t) . 

The function s w (t) is given below for each a. 

3 ALGORITHMS FOR EDF CALCULATION 

Our purpose is to obtain practical numerical approximations of (6). We give two versions of the 
algorithm: the simplified version merely truncates the sum in the exact formula; the full version 
maintains the number of summation terms below a presassigned threshold and avoids catastrophic 
roundoff errors. They have the same inputs, output, function definitions, and initial step. Some 
explanation of the approximations is given in the Appendix. Because the results are invariant to 
time scaling, we may set r = 1, To = 1/m. 

All arithmetic is to be carried out in double- precision floating point. Operations that give signed 
integers are the floor function |_xj (greatest integer that is < x) and ceiling function \x~\ = — [— x\ 
(least integer that is > x). 
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3.1 Inputs 

a = frequency noise exponent 
a = 2, 1,0, -1,-2, -3,-4 

Noise type = WHPM, FLPM, WHFM, FLFM, RWFM, FWFM, RRFM 
d = order of phase difference 

d = 1: first-difference variance (included for completeness) 
d=2: Allanvariance 
d = 3: Hadamard variance 
Restriction: a + 2d > 1 
to = averaging factor t/tq, positive integer 
F = filter factor 

F = 1: modified variance 
F = to: unmodified variance 
S = stride factor (estimator stride = t/S) 

S = 1: nonoverlapped estimator 
S = to: overlapped estimator 
N = number of phase data with sample period tq 


3.2 Output 

edf = equivalent degrees of freedom of the variance estimator 


3.3 Constant and function definitions 

Set an integer constant J max (used only in the full version); we suggest J max = 100. 

The formal arguments of the following functions have the same names as the input arguments of 
the main algorithm. 


1. Define the function s w (f, a) as follows: 

a 2 1 

s w (t,a ) -\t\ t 2 In \t\ 


-1 

— t 4 In I 


The entries with In If I must evaluate to 0 when t = 0. 


-2 -3 -4 


|5 +6 


t b In | 


2. Define the function 


s x (f, F, a) = F 2 A 1 /f A_ 1/f s w (f, a) 


= F 1 


1 


2 Syj (t, Q() Syj ( t ^ , QL I S W ( t + 


1 


s x (f, oo, a) = s w (t, a + 2) , —4 < a < 0. 

3. Define the function 

s- (f, F, a, d) = (AiA_i) d s x (f, F,a), d= 1, 2, 3; 


(7) 


( 8 ) 


(9) 
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that is (with dependence on F and a suppressed on the right sides), 

s* (t, F, a, 1) = 2 s x ( t ) - s x (t - 1) - s x (t + 1) , 

s z (t, F, a, 2) = 6s x ( t ) - 4 s x (t - 1) - 4s* (f + 1) + s* (t - 2) + s x (t + 2) 

s z (t, F, a, 3) = 20s* (t) - 15s* (t - 1) - 15s* (t + 1) 

+ 6s* (t — 2) + 6s* (f + 2) — s* (t — 3) — s* (i + 3) . 

4. Define the function 


BasicSum ( J, M, S, F. a, d) = s 2 z (0, F, a, d) + ^1 — — \ s 2 z ( — ,F,a,d 


+ 2 ^( v 1 MJ s z[ s ,F,a,d 

3 = 1 


3.4 Initial steps for both versions 

1. Compute M, the number of summands in the estimator, as follows: 

L = — + md (an integer), 

F 


M = 1 + 


S (N — L) 


if N > L, otherwise there are not enough data. Remark: L is the length of the filter applied to the 
phase samples. 


2. Let 


J = min (M, (d +1)5). 


3.5 Main procedure, simplified version 


This is just one step: 


edf s 2 z (0, F, a, d) M 


BasicSum ( J, M, S,F, a, d) . 


To check the effect of the truncation, one can also try a larger value of J, say min (M, 65). 


3.6 Main procedure, full version 

T M 

Let r = 

There are four cases. The calculations use coefficients from three numerical tables. 
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3.6.1 Case 1. Modified variances: F = 1, all a 


This case also applies to unmodified variances when F = m = 1. 

d I I '.WAX 

!, = 2 77, -I 1 — w BasicSum (J,M,S,l,a, d) 

edl (0, l,a,d)M 

Else if J > J max and r > d + 1, take ao,oi from Table 1; then 

1 1 / ai \ 


= ( «o 


edf r V 

Else let m! = max (not necessarily an integer); then 
r 


edf s\ (0, 1, a, d) </„ 


-BasicSum (J max , J max , m , 1, a, d) 


3.6.2 Case 2. Unmodified variances, WHFM to RRFM: F = m, a < 0 

If 'I 'i AIAX 

If m (d + 1) < J max then let m! = m else let m! = oo. Then 

iS = ,2(0. m'! n ,rf)M BaSicSum (• J ' S - m< ’ “■ 

Else if J > J max and r > d + 1, take ao,oi from Table 2; then 

1 1 / Ol \ 


edf r 


= - o 0 


Else let m' = (not necessarily an integer); then 

r 


edf s 2 z (0, oo, a, d ) J n 


-BasicSum ( J max , J max , m\ oo, a, d) 


3.6.3 Case 3. Unmodified variances, FLPM: F = m, a = 1 

d '1 d max 


edf (0,m,l,d)M 


BasicSum ( J, M, 5, m, 1, d) . 


Remark: For this case, m must be less than about 10 6 to avoid roundoff error. 

Else if J > J max and r > d + 1, take ao,ai from Table 2 (cc = 1), bo,bi from Table 3; then 


edf (bo + bihun) 2 r 


Else let m' = (not necessarily an integer); then 

r 


edf (bo + b\ In m) 2 J n 


-BasicSum ( J max , J max , m! , m', 1, d) 
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3.6.4 Case 4. Unmodified variances, WHPM: F = m, a = 2 

This calculation is exact 

binomial coefficient — — — 7 . 

k\ (n — k)\ 

Let K = [?’] . 

If K < d 


, and can be expressed in closed form. In these formulas, 


n\ 


Else 


1 

edf 


1 

M 


1 + 



K—l 


E 1 


k= 1 


k\ f 2d \ 2 

rj V d-kj 


11/ _ ai 

edf M V ° r 


where 


cio = 




also given in Table 2 (a = 2). 


3.7 Tables 


Table 1. Coefficients for modified variances. 



d 

= 1 

d 

= 2 

d 

= 3 

a 

ao 

a± 

ao 

«i 

a o 

ai 

2 

1 

2 

3 

0.840 

l 

3 

0.345 

7 

9 

0.997 

1 

2 

0.616 

22 

25 

1.141 

‘2 

3 

0.843 

0 

1.079 

0.368 

1.033 

0.607 

1.184 

0.848 

-1 



1.048 

0.534 

1.180 

0.816 

-2 



1.302 

0.535 

1.175 

0.777 

-3 





1.194 

0.703 

-4 





1.489 

0.702 


denotes the 
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Table 2. Coefficients for unmodified variances. 



d 

= 1 

d 

= 2 

d 

= 3 

a 

(io 

fll 

do 

a i 

do 

ai 

0 

3 

1 

35 

1 

231 

3 

Zj 

2 

2 

18 

100 

2 

i 

78.6 

25.2 

790. 

410. 

9950. 

6520. 

0 

2 

l 

2 

l 

7 

l 

3 

6 

3 

3 

9 

2 

-1 



0.852 

0.375 

0.997 

0.617 

-2 



1.079 

0.368 

1.033 

0.607 

-3 





1.053 

0.553 

-4 





1.302 

0.535 


Table 3. Coefficients for logarithmic denominator, unmodified variances, FLFM (a = 1). 

d = 1 d = 2 d = 3 
fro h b 0 b\ b 0 bi 

6 4 15.23 12 47.8 40 


4 EXAMPLES 

First, we must point out that the new edf values differ somewhat from older ones because of our 
choice of phase noise model. Previously (for a < — 1) the phase was generally assumed to have 
a pure f a ~ 2 spectrum; here, the phase is modeled as the first difference of an f a ~ 4 process. For 
example, consider overlapped Allan variance, white FM, 1025 phase data. The old results are from 
[ 14 ] (close to those of [ 7 ] . 


t/t 0 

1 

2 

4 

8 

16 

32 

64 

128 

256 

512 

old edf 

682 

584 

354 

186.3 

93.4 

45.8 

21.8 

9.83 

4.01 

1 

new edf 

801 

554 

314 

170.0 

88.5 

44.4 

21.8 

9.83 

4.00 

1 


The old and new results are in practical agreement at the larger values of r, where the results 
matter more. By allowing this mild discrepancy, we were able to make the algorithm simpler and 
more uniform. 

Figures 1 and 2 (at the end of the paper) show examples of edf, computed by the new algorithm, 
as a function of noise type and of variance type with other parameters fixed. 

5 CONCLUSION 

The stability variances considered here can all be put into the same form, namely, the mean- 
square average of the output of a finite-difference filter acting, not on the phase samples, but on 
their cumulative sums. With this insight, we have been able to make an algorithm for computing 
the equivalent degrees of freedom of variance measurements. It covers all the commonly used 
variances and estimators, and then some. There is now a single unified source for these uncertainty 
calculations instead of the many papers that each cover only one variance and perhaps only one 
estimator of it. Complete pseudocode for the new algorithm is given here; software is available on 
request [ 23 ]. 
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6 APPENDIX: EXPLANATIONS OF APPROXIMATIONS 

As is, (6) is unfit for numerical computation. We find empirically that s\ (t) tends rapidly to zero 
as t increases beyond d. For the accuracy needed here (a few percent), there is no point in allowing 
j/S to be more than d+ 1. Indeed, for sufficiently large t the calculation of s\ (t) blows up from 
roundoff error, even in double precision, because linear combinations of large s w values are taken 
to get small s z values. At the very least, one should truncate the sum at j = (d + 1) S , as in the 
simplified version of the algorithm. 

The full version of the algorithm uses the following general strategy. If J < J max we do the 
summation (6). When J > J max there are two cases. First, if M > (d+l)S then S = m > 
Jmax/ ( d + 1) >> 1. We truncate the sum at (d + 1) S and approximate it by an integral; this gives 


where 


1 

edf Vd 



( X “ /) s z (*) dt 




M 

r = ~s 1 a ° 



dt, 


a i 



s\ ( t ) tdt. 


These coefficients can be evaluated in advance. Second, if M < (d + 1) S then we do another 
summation in which J is reduced from M to J max and S is reduced proportionately from m. 


The extra term for j = J in BasicSum makes the sum a trapezoidal approximation to the integral, 
whether or not the sum is truncated. 


This method works straightforwardly for Case 1; indeed, in this setting the modified variances are 
simpler than the unmodified ones. In Case 2, when m is large we compute s z (t) using the limiting 
form s x (t, oo), which is actually —s'/, (t). This means that we are treating x (t) as w' (t), the process 
w (t) being differentiable in the mean-square sense. 

The most troublesome case is the overlapped estimators of the unmodified variances for flicker PM. 
As S = m — > oo, s z (t) approaches a function with logarithmic singularities. The factor bo + b\ lnm 
is an asymptotic form of s z (0). It would be possible (though inconvenient) to add another large-m 
subcase as in Case 2, but one does not expect flicker PM to be the dominant noise type when m is 
large. 

Case 4 is constructed by knowing that the phase samples are accurately uncorrelated when w (t) 
is a Wiener process. The simplified computation (13) is correct, but wasteful, because s z (t) is a 
linear combination of hat-shaped peaks of width 2/m. 
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Overlapping ADEV EDF for N=100 



Figure 1. Edf vs. averaging factor with power-law noise type as a parameter: Allan variance, 
overlapped estimator, 100 phase samples. 



1 10 100 
Averaging Factor, m 


Figure 2. Edf vs. averaging factor for three stability variances: overlapped estimator, white FM, 
100 phase samples. 
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Questions and Answers 

DON PERCIVAL (University of Washington): Your algorithm, your recipe handles the integer power 
laws. How hard would it be to include something like the power law of minus five-thirds, which would make 
people in the atmospheric turbulence community very happy. 

CHUCK GREENHALL: It would not be difficult to do a specific extra power. It would be more work to 
include a whole continuous range. The algorithm is split up into four cases for purposes of numerical 
computation. So it would be a modest computation. 

MARC WEISS (National Institute of Standards and Technology): I guess I missed it, but it sounds like 
you need to know the power law before you can get the confidence. 

GREEN HALL : That is correct. 

WEISS: So how do you suggest doing that? 

GREEN HALL : Well, a couple years ago there was a paper in which both of us were co-authors, and it 
actually published the method being in the Stable software. Bill and I have a paper coming up at EFTF in 
which we have developed a simpler and better method involving the lag-1 autocorrelation. 

DEMETRIOS MATSAKIS (U.S. Naval Observatory): Maybe a better way to phrase the other question 
was: what percentage of the error is due to the error bars? You have a certain contribution because you do 
not know how to fit your power law. So how much does it raise your error? 

GREENHALL: Well, I don’t know. We were pulling ourselves up by our bootstraps here, and you’re 
probably most conservative to assume the more red power law if you are uncertain about it. That is the best 
thing I can say. 
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